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Abstract: The matrix model formulation of superstring theory offers the possibility to 
understand the appearance of 4d space-time from lOd as a consequence of spontaneous 
breaking of the SO(IO) symmetry. Monte Carlo studies of this issue is technically difficult 
due to the so-called sign problem. We present a practical solution to this problem gen- 
eralizing the factorization method proposed originally by two of the authors (K.N. A. and 
J.N.). Explicit Monte Carlo calculations and large- extrapolations are performed in a 
simpler matrix model with similar properties, and reproduce quantitative results obtained 
previously by the Gaussian expansion method. Our results also confirm that the sponta- 
neous symmetry breaking indeed occurs due to the phase of the fermion determinant, which 
vanishes for collapsed configurations. We clarify various generic features of this approach, 
which would be useful in applying it to other statistical systems with the sign problem. 
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1. Introduction 

One of the biggest puzzles in string theory is that typical space-time dimensionality turns 
out to be higher than the macroscopically observed four dimensions. As a possible approach 

to this problem, one may think of space-time as an emergent notion, which arises effectively 
from matrix degrees of freedom. This idea of "emergent space-time" is nicely realized in the 
gauge-gravity duality including the famous example of the AdS/CFT correspondence [1]. 
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Recently the gauge-gravity duality has been demonstrated from first principles by 
Monte Carlo simulations of super symmetric gauge theory in one dimension with 16 super- 
charges [2-9]. (See ref. [10] for a review.) This theory can be obtained formally by the 
one-dimensional reduction of ten-dimensional J\f = I super Yang- Mills theory, and it 

provides a low energy description of a stack of N DO branes in type IIA superstring theory. 
In particular, in the iV — t- oo limit with large 't Hooft coupling, the one-dimensional gauge 
theory was conjectured to be dual to the 0-brane solution in type IIA supergravity [11]. 
Monte Carlo calculations have been performed on the gauge theory side for various quan- 
tities such as the internal energy [4-7], the Wilson loop [8] and the two-point correlation 
functions [9], and the results provided first-principle confirmation of corresponding pre- 
dictions based on the gauge-gravity duality. One of the most important conclusions [6] is 
that the black hole thermodynamics of the 0-brane solution has been understood micro- 
scopically in terms of open string degrees of freedom attached to the DO branes, which are 
described by the gauge theory. In other words the Id U(A^) supersymmetric gauge theory 
describes the interior structure of the black hole in ten dimensions. 

As a closely related but more ambitious conjecture, the IIB matrix model was proposed 
as a nonperturbative definition of type IIB superstring theory in ten dimensions [12] . This 
model can be formally obtained by the zero-dimensional reduction of ten-dimensional = 
1 U(A^) super Yang-Mills theory, and it has manifest SO(IO) symmetry. Here the space-time 
is represented by the eigenvalue distribution of the 10 bosonic matrices [13]. Therefore, 
the model offers the possibility to realize dynamical compactification,^ where the extra 
dimensions become small due to the spontaneous symmetry breaking (SSB) of SO(IO). 
This scenario has been supported by explicit calculations based on the Gaussian expansion 
method (GEM). By comparing the free energy of the SO{d) symmetric vacua with d = 
2, 4, 6, 7, it was shown that d = 4 is favored at the 3rd order of the expansion [15]. Higher 
order calculations confirmed this conclusion [16]. (Sec rcfs. [17,18], however.) 

Just as Monte Carlo simulations have been useful in studying the gauge-gravity du- 
ality from first principles, they are expected to be useful also in addressing the issue of 
dynamical compactification in the IIB matrix model and many others. Indeed early works 
on zero-dimensional matrix models [19-23] have provided us with a wealth of information 
on the large- limit of the models and the nonperturbative dynamics of their degrees of 
freedom. When one applies such an approach to the IIB matrix model, one encounters 
a serious technical problem called the sign problem since the determinant (or PfafBan, 
strictly speaking) one obtains from integrating out the fermionic matrices is complex. Due 
to the fluctuation of the phase, huge cancellations occur in Monte Carlo integration over the 
bosonic matrices. To overcome this problem, a promising method termed the factorization 
method has been proposed [24], tested on simple models [25,26] and used in simulations 
of finite density QCD [27]. (See refs. [28,29] for related works on the QCD phase diagram 
and refs. [30] for other approaches to the complex action problem.) It resembles previously 
proposed density of states methods [31] but suggests more general and efflcient ways of 
slicing the configuration space. With the proposed technique it is possible to compute the 

^See ref. [14] for an earlier work suggesting this possibility in the framework of bosonic string field theory. 
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chosen distribution function in the region which is enhanced due to relatively small fluc- 
tuations of the phase. Expectation values can be computed as the easily located minima 
of the free energy thereby drastically reducing the propagation of errors. It also has the 
merit of being in principle applicable to a wide range of models and its possible success 
can have a great impact on various physical problems. 

In fact the fluctuations of the phase of the fermion determinant is expected to play a 
crucial role in the dynamical conipactification in the IIB matrix model [32].^ It was shown 
in ref. [32] that the phase vanishes for collapsed configurations, and that the enhancement 
due to this property can compensate the entropic suppression for such configurations. In 
ref. [24] it was shown that the competition between the two effects can make the dom- 
inant configurations have very different length scales than in the phase-quenched model 
suggesting a mechanism for dynamically generating small and large dimensions. Models 
without the phase factor have also been studied by Monte Carlo simulations, which pro- 
vided strong evidence that the SSB does not occur in such models [20-22]. These results 
stress the importance of the role played by the phase factor. 

In this paper we study a simple zero-dimensional matrix model [34], which realizes 
the dynamical compactification. The (non-supersymmetric) model has SO (4) rotational 
symmetry and consists of four N x N bosonic matrices and A'^^ flavors of Weyl fermions in 
the fundamental representation of SU(A^). The integration over fermions yields a complex 
determinant. In ref. [34], the large- A'^ limit is taken by keeping the ratio r = Nf/N fixed, 
and it is shown analytically for infinitesimal r that the SO (4) symmetry is broken down 
to S0(3). If the phase is quenched, the SSB is shown not to occur. For finite r, the 
model has been studied using the GEM [35] . While the exact results for infinitesimal r are 
consistently reproduced, it is shown that the symmetry actually breaks down to S0(2) at 
finite r due to the stronger effect of the phase. 

Motivated by these results, we perform Monte Carlo simulation of this simplified model, 
which is expected to serve as a testing ground for the ideas discussed above concerning the 
SSB in the IIB matrix model as well as for the method of simulating systems with a com- 
plex action. Besides being lower dimensional, the model is much easier to simulate than the 
IIB matrix model with the computational effort increasing as instead of AT^, since the 
fermions are in the fundamental representation instead of the adjoint representation. We 
find that the method with appropriate generalization samples efficiently the configuration 
space, heavily suppressed in the phase-quenched model but important for the full model, 
overcoming the so-called overlap problem. The distribution function and the phase of the 
fermion determinant exhibit nice scaling properties which are important for the extrapo- 
lation to large N and to the regions of configuration space in which the phase fiuctuation 
obscures the Monte Carlo measurements. We present detailed analyses of this scaling, 
and provide its clear understanding based on simple theoretical arguments. Thus we have 
confirmed that the SSB indeed occurs in this model due to the effect of the phase. 

We first apply the factorization method to a single observable as originally formulated 
in ref. [24]. While the results show that the method works to some extent, we also realize 

^See refs. [13,33] for discussions on other possible mechanisms. 
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certain discrepancies from the GEM results as well as some puzzles, which we attribute to 
the remaining overlap problem. In order to solve this problem, we generalize the method 
to multiple observables, and show that the GEM results can be consistently reproduced. 
The importance of controlling multiple observables in the factorization method is reported 
briefly in our previous publication [36]. 

This paper is organized as follows. In section 2 we review the basic properties of the 
model, and discuss the sign problem that arises in its Monte Carlo studies. In section 3 
we investigate the model by the factorization method with a single observable. In section 
4 we generalize the factorization method to multiple observables. Section 5 is devoted to 
a summary and discussions. In appendix A we describe the details of the algorithm used 
in our Monte Garlo simulation. In appendix B we discuss some properties of the fimctions 
that play an important role in the method. In appendices G and D we present the details 
of the large- extrapolations made in our analyses. 

2. The model and the sign problem 

The model we study in this paper is defined by the partition function [34] 

Z = J dAdiljdtpe-^^''+^'\ (2.1) 

5b = ^iVtr(^^)2 , (2.2) 

Sf = -iPiir^)apA^i;f^. (2.3) 

The bosonic degrees of freedom are represented hy N x N Hermitian^ matrices (/x = 
1, • • • ,D). Hereafter we assume D to be even.'^ The fermionic degrees of freedom are 

represented by -{[ja and ipa, which have a spinor index q = 1, • • • ,p, where p represents the 
number of components of a Z)-dimensional Weyl spinor 

p = 2^/2-1 . (2.4) 

They also have a flavor index / = 1, • • • ,iVf, where iVf represents the number of flavors. 

— J" J- 

The fermionic variables ipa and ipa are iV-dimensional row and column vectors, respectively, 
so that the actions (2.2) and (2.3) have an SU(A^) symmetry. The p x p matrices are 
SO{D) gamma matrices after the Weyl projection. Thus the actions (2.2) and (2.3) have an 
SO{D) symmetry, where the bosonic variables A^ transform as a vector and the fermionic 
variables transform as Weyl spinors. The fermionic part can be regarded as the zero volume 
limit of Weyl fermions interacting with a background gauge field via fundamental coupling. 
Integrating out the fermions, we obtain 

Z = j dAe-^^Zi[A] , (2.5) 



^In rcf. [35] the matrices are assumed to be also traceless to simplify the calculations in the Gaussian 
expansion. In our Monte Carlo studies, we do not impose the tracelessness condition to avoid unnecessary 
complication, since it is irrelevant in the large- AT limit. 

*For odd D, the model can be defined using Dirac fermions instead [34]. In that case, however, the 
fermion determinant is real, and no SSB of SO(Z)) is expected. 
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where Zf[A] = (detD)^f and V = T^A^ is a pN x pN matrix. Let us then discuss the 
properties of the fermion determinant det T> for a single flavor. First of all, it is complex 
in general. Under parity transformation Ad ^ —Ad, Ai Ai [i ^ D), it transforms as 
detP (detP)*. This implies that detX> is real^ for configurations with Ad = 0. Prom 
this fact alone, it follows that the phase of the determinant becomes stationary for config- 
urations with Ad = Ad-i = since one cannot have a phase fluctuation within a linear 
perturbation around such configurations.® These properties of the fermion determinant are 
analogous to those found in the IIB matrix model [32]. 

As in the case of the IIB matrix model [32] , one can extend the above argument further 
to arrive at the following statements, which actually play an important role in our analysis. 
Let us define "d-dimensional configurations" {d > 1) as such configurations A^ that can 
be transformed into ^d+i = ■■■ = Ad = hy sm SO{D) transformation. Then for general 
d- dimensional configurations, we have 

^'^Q^ =0 for k = l,--- ,{D-l)-d, (2.6) 

where P represents the phase of the fermion determinant Zf [A] . 

We take the large- iV limit with r = N^/N fixed, which corresponds to the Veneziano 
limit. ^ This is needed to make the fermionic degrees of freedom contribute to the partition 
function comparably to the bosonic degrees of freedom in the large-iV limit. Whether the 
SSB of SO(L') occurs in that limit is the issue we would like to address. Por that purpose, 
we consider the "moment of inertia tensor" [13,20] 

T^, = ^tT{A^A,) , (2.7) 
and its real positive eigenvalues (n = 1, • • • , D) ordered as 

Ai > A2 > ••• > Ad . (2.8) 

The vacuum expectation values (VEVs) of these eigenvalues (A„) play the role of the order 
parameters. If they turn out to be unequal in the large-A^ limit, it implies the SSB of 
SO{D). In this model, the sum of the VEV of all the eigenvalues is given exactly as^ 

Y.{\n) = (^^tr{A,f\=D+pT . (2.9) 



^In this case, the fermion determinant Zi[A\ for A'f flavors is actuaUy real positive for even A^f, although 
it is not necessarily so for odd iVf. There is some numerical evidence, however, that configurations with 
positive determinant dominate statistically in the large- iV limit. Therefore we consider that there is no 
distinction between odd N{ and even N{ in the large-AT limit. In the present work, we always use even Nf 
to avoid this subtlety. 

^Obviously the same statements hold for configurations obtained by S0(i3) rotations. 

7pQj. ^ _ ^]^g fermionic variables can be written in terms of N x N matrices ('l'a)i/ and {'^a)fi- 
The fermionic part of the action becomes St = —{T^)aptr{^aA^'i p), which can be compared with a term 
S{ = -{ri^)aptT{^c,[Aft, *;3]) in the IIB matrix model. 

®This can be derived by rescaling A'^ — k.2 A^. Defining Z(k) — f dAZ{[A] exp [— iiVtr(^Jj)^] and using 
Zi[A'] = KP'-^'/2^f[A], we obtain = Finally {j^tr{A^f) = ^^^logZ(fi)\ 

I K — 1 

D + pr. The same derivation with the replacement Z{[A] i-> |2f holds for the phase-quenched model. 



-5- 



From now on, let us consider the D = 4 case, which imphes p = 2 due to (2.4). The 
gamma matrices are given by 



Ti = (71 = 




, r2 = (72 




Using a simple technique familiar in random matrix theory, it was shown in the large- 
limit and at infinitesimal r that the VEV (An) is given as [34] 



(An) = 



(2.10) 



1 + r + o(r) for n = 1, 2, 3 , 
1 — r + o(r) for n = 4 , 

which implies that the S0(4) symmetry is broken down to S0(3). This SSB is associated 
with the formation of a condensate (V'aV'a)) which is invariant under S0(3) only. 
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Figure 1: The VEVs (A„)o (n = 1, 2, 3 and 4) in the phase-quenched model (2.11) are plotted for 

r = 1 (Left) and r = 2 (Right) against 1/A^. The data for iV > 8 can be nicely fitted to straight 
lines meeting at the same point (1 + r/2) at iV = oo, which demonstrates the absence of SSB for 
each r. 

Here the effect of the phase plays a crucial role. To see that, let us consider the 
"phase-quenched model" 

Zo = J ci^e-^°[^l , (2.11) 

So[A] = S^[A] - Nf log I detV[A] \ . (2.12) 

Eq. (2.9) holds for this model as well (See footnote 8.), and therefore the absence of SSB 
would imply in the large- A?^ limit 

(An)o = l + ^ for all n = 1,2,3,4, (2.13) 

where the VEV ( • )o is taken with respect to (2.11). Indeed this can be confirmed at 
infinitesimal r [34]. Here we simulate the system (2.11) for r = 1, 2 and plot the eigenvalues 
(A„)o against 1/A^ in fig. 1. We find that the data can be nicely fitted to straight lines, 
and the extrapolation to = oo gives the value 1 -|- r/2 for all n = 1, 2, 3 and 4 in accord 
with (2.13)9 with accuracy better than 0.3%. This implies, in particular, that the SSB of 
^Similar results have been obtained for r = 0.25, 0.5, 1.5, 4.0. 
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so (4) does not occur in the phase-quenched model. 

In what follows we will study the normalized eigenvalues 

An 



'^n — 



(A. 



n/O 



(2.14) 



The deviation of (A„) from 1 represents the effect of the phase. The relevant question is 
whether this deviation depends on n in the large- iV limit. 

At finite r < 2, the large- iV limit 
of the full model (2.1) was studied by 
the GEM up to the 9-th order [35] in 
the same way as it was used in the JIB 
matrix model [15, 16]. The results for 
r = 1 and r = 2, which are the cases 
we focus on in this paper, are summa- 
rized in Table 1. We have translated the 
results for the unnormalized quantities 
(A„) in ref. [35] to the normalized ones 
(An) using (2.13) and (2.14). 
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1.17 


1.4 


1.23 


1.7 


(A2) 


1.17 


1.4 


1.23 


1.7 


(A3) 


1.17 


0.7 


1.23 


0.5 


(A4) 


0.5 


0.5 


0.31 


0.1 


free energy 


-1.5 


-1.8 


-1.9 


-3.6 



Table 1: The results for the normalized eigenvalues 
(A„) and the free energy obtained by the GEM at 
N = 00 with the S0(3) and S0(2) ansatz [35]. 



The free energy obtained with the S0(2) ansatz was found to be smaller than the one 
with the S0(3) ansatz, which implies that the true vacuum is only S0(2) symmetric. At 
smaller r [35], the free energy and the VEV (An) obtained with the S0(2) ansatz asymptotes 
to the ones obtained with the SO (3) ansatz, which ensures the consistency with the analytic 
results (2.10) obtained at infinitesimal r. 

Monte Carlo simulation of the full model is not straightforward since the complex 
integrand of (2.5) cannot be regarded as the Boltzmann weight, and hence the idea of 
importance sampling is not applicable. One approach is to rewrite (2.5) as 



and to calculate expectation values by reweighting 

(Oe^^)o 



{O) 



(2.15) 



(2.16) 



where the VEVs on the right-hand side are taken with respect to the phase-quenched model 
(2.11), and hence can be evaluated, in principle, by standard Monte Carlo techniques. The 

VEV (e*'")o is nothing but the ratio of the partition functions Z/Zq, and it decreases 
exponentially at large as e~^^^^, where AF > is the difference in free energy of the 
two systems. This can happen due to huge cancellations from the factor e*^. As a result, 
one needs 0(e'^°"®*-^^) configurations to compute an observable with given accuracy by 
using the formula (2.16) directly. This is called the sign problem. 

Another closely related problem is the so-called overlap problem. By using Zq in 
sampling configurations, we usually sample the wrong region of configuration space, which 
has little overlap with the dominant configurations in the system given by Z. This problem 
becomes exponentially hard at large N and makes importance sampling inefficient. 
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3. Factorization method with a single observable 

In this section we perform Monte Carlo studies of the model (2.5) by applying the factor- 
ization method to a single observable as originally proposed in rcf . [24] . While the analysis 
shows various encouraging behaviors, we will see that the method is only partially suc- 
cessful in reproducing the known results of the GEM. We argue that this is due to the 
remaining overlap problem, and discuss in section 4 how it can be solved by generalizing 
the method to multiple observables. 

3.1 the basic idea 

Let us consider calculating the VEV (A„) by Monte Carlo simulation. Instead of using the 
reweighting formula (2.16) directly, we first rewrite the VEV as 

/•oo 

(An) = / dxxpn{x) (3.1) 
JO 

in terms of the distribution function 

Pn{x) = (5{x-\n)) . (3.2) 

Applying now the reweighting (2.16) to the right-hand side, one finds that it factorizes as 

Pn{x) = ^P^^\x)Wn{x) . (3.3) 

The real parameter C is a normalization constant given by^'^ 

C-M {e'^)^ = {cosV)o , (3.4) 

which need not be calculated in the present method. The real function pn^ {x) is nothing 
but the distribution function of A„ in the phase-quenched model defined by 

p^^\x)^^{5{x-K))^, (3.5) 

which is peaked at a; = 1 due to the chosen normalization (2.14). The function Wn{x) in 
(3.3) is defined by 

Wn{x) =^ {e''^)n,x = {cOsT)n,x , (3.6) 

where ( • )n,x denotes a VEV with respect to the partition function 

Zn,x = j dAe-^o 5{x - \n) . (3.7) 



It turns out that Wn{x) > 0, which simplifies our analysis significantly. 



11 



^''In the second equality, we have used the fact that the phase F flips its sign under the parity transfor- 
mation Ad i— >■ —Ad, which is a symmetry of the phase-quenched model (2.12). A similar comment applies 
also to the second equality of (3.6). 

^^This property does not always hold in applications of the factorization method to a system with the 
complex action problem. See ref. [25] for a study in such a case. 
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Using the saddle point approximation in the integral (3.1), the problem of determining 
(A„) can be reduced to that of minimizing the "free energy" 

J^n{x) = - \ogpn{x) , (3.8) 

which simply amounts to solving 

fi'\x) Aiogp(0)(^) = -A.iog^^(^) . (3.9) 

It turns out that both sides of this equation scales as 0{N^). Therefore, the saddle point 
approximation actually becomes exact in the N oo limit. 

3.2 practical implementation of the method 

As we have discussed above, the calculation of the VEV (A„) reduces to the determination 
of the two functions u;„(x) and fn\x). Here we discuss how we can measure them by 
Monte Carlo simulation following ref. [24]. For that purpose we approximate the delta 
function S{x — \n) in (3.7) by the Gaussian function and simulate the system defined by 

Zn,v = J d^e-^^°+^(^")> , V{z) = \liz- 0' , (3.10) 

where 7 and ^ are real parameters. Let us then consider the distribution function of A„ for 
the system (3.10), which is given by 

Pn,v{x) = {6{x - An))^^ OC p^^\x) exp { - v[x{Xn)o)} , (3-11) 

where ( • )n,v represents a VEV with respect to Zny. The position of the peak, which we 
denote by x-pj Cciri be obtained by solving 

= ^log pn,v{x) = /^°)(x) - {\n)oV'{x{Xn)o) • (3.12) 

Since the distribution function pny{x) is sharply peaked at Xp for sufficiently large 7, 
we can use the VEV of A„ as an estimator of Xp, i.e., 

Xp = {Xn)n,V ■ (3.13) 

By varying the value of ^ in (3.10), we obtain the functions Wn{x) and fn\x) by^^ 

WniXp) = {cOSr)n,V , (3.14) 
/W(Xp) = (A„)o V'(^{Xn)n,v) = 7 (An)0 {{K)n,V - ■ (3-15) 

The parameter 7 should be chosen large enough to make the fluctuation of A„ smaller than 
the required resolution in x. For the purpose of evaluating fn\x), however, one should not 



^^One might naively tiiink tiiat the function pn\x), and hence the function fn \x), can be obtained 
by measuring the distribution of A„ in the phase-quenched model Zq. By such a direct method, however, 
accurate calculation is possible only in the vicinity of the peak x = 1. Note that the solution of (3.9) is 
typically not very close to the peak. 
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use too large 7. Theoretically, the right-hand side of (3.15) should converge to the correct 
value in the 7—7-00 limit, which implies that (^{Xn)n,v ~ I/7. Hence a small error in 

{'^n)n,v propagates to fn\x) by the factor of 7. 

Note that the computation of Wn{x) based on (3.14) is hard due to the sign problem. 
However, we only have to obtain Wn{x) in the region where the (possibly local) minima of 
the free energy (3.8) are likely to exist. This typically implies the region where Wn{x) is 
not extremely small. In such a region, one can make a sensible large- AT extrapolation 

^n{x) = lim -^logw;„(a;) . (3.16) 

AT— >-oo iV 



Then one can use the equation 



= -^^n(x) (3.17) 



to obtain the VEV (A„) . Since there is no sign problem in the calculation of fn^ (x) , one 
may obtain the left-hand side of (3.17) at larger (if necessary) than the values of N 
used to obtain $„(x) in eq. (3.16). Note also that the error on both sides of (3.17) due to 
statistics and finite N does not propagate exponentially to (A„) as a direct computation 
based on (3.1) would imply. 

3.3 asymptotic behaviors of the functions Wn{x) 

In our approach the effects of the phase T are represented by the right-hand side of (3.17). 
Therefore, the large- A?^ extrapolation (3.16) is the most crucial step for the success of the 
method. In the left columns of figs. 2 and 3, we plot our results for log Wn{x) for r = 1 
and r = 2, respectively. 

As one can see from these plots, the function Wn{x) approaches 1 for x <^ 1 and/or 
x ^ 1. These are the regions, which are favored by the effects of the phase. We can un- 
derstand these properties as follows. Let us recall that Wnix) defined by (3.6) is the VEV 
of e^^ in the ensemble (3.7). The dominant configurations in the ensemble for x <C 1 have 
(5 — n) shrunken directions due to the ordering (2.8), and hence they are approximately 
configurations with dimension d = {n — 1). Similarly, dominant configurations for x ^ 1 
have n extended directions, and hence they are approximately configurations with dimen- 
sion d = n. Since the phase of the determinant vanishes for collapsed configurations (i.e., 
with 1 < d <3 dimensions) as we explained below eq. (2.5), the VEV of e*^ approaches 1 
when such configurations dominate in the corresponding ensemble. 

Moreover, we can deduce the asymptotic behaviors of the functions Wn{x) as 

1, .^ J-Cnx^- (a;« l,n = 2,3,4) , 
^log^.(-) - I _^„,-(4-„) ^ 1,2,3) . ^^-''^ 

This can be derived from the property (2.6) of the phase. It follows that the fluctuation 
of the phase is of the order of 6T {SA/\A\)^^'^ around a collapsed configuration with 
1 < d < 3 dimensions, where SA and \A\ represent the typical scale of A^ in the shrunken 
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and extended directions, respectively. Due to the definition (2.7), it is expected that 5y4/|^| 
is proportional to \fx (and 1/v^) for a; <C 1 (and a; 3> 1), respectively. When a; <C 1, the 
distribution of the phase is therefore expected to have a width a oc (y^)^"". Assuming 
that the distribution is Gaussian/^ we can evaluate the expectation value of e^^ by using 
the formula [28] 



/ 



Thus we obtain — logu'„(a;) = ^a^ oc x^^". Similarly, when x 3> 1, the distribution is 
expected to have a width of the order of (l/^/x)^^"', and hence — logWnix) oc x~^'*'~'^\ 

We fit our data to the asymptotic behaviors (3.18) and extract the coefficients and 
dn as described in appendix C. It turns out that the finite-A?^ effects in these coefficients 
are of the order of 1/N. Based on this observation, we make large- A" extrapolations to 
obtain (3.16) in the asymptotic regions, which is also plotted in the left columns of figs. 2 
and 3. The two solid lines represent the margin of error due to the uncertainty from the 
large- AT extrapolation. 

The regions of x in which ^n{x) becomes small represent the regions in which the 
sign problem becomes severe. As we approach such regions, it becomes more difficTilt to 
obtain data points, in particular, for larger A^. A big virtue of the factorization method is 
that we can use the function <&n(a^) in the asymptotic regions to extrapolate towards the 
region where the data points are not available due to the sign problem. One should keep 
in mind, however, that the true function ^n{x) may deviate from the asymptotic behavior 
as one approaches x ~ 1, which causes certain systematic error. We will discuss this issue 
in sections 4.4 and 4.6. 

3.4 results for (A^t) 

On the right columns of figs. 2 and 3, we plot our results for the function ■^fn\x) for 
r = 1 and r = 2, respectively. Since there is no sign problem in this calculation, we can 
obtain results for much larger N, which clearly show the large-AT scaling behavior. Note 
also that the function ■^fn\x) crosses zero at x ~ 1 corresponding to the fact that Pn\x) 
is peaked at x ~ 1. The asymptotic behaviors of the functions fn\x) in the regions x ^ 1 
and X ^ 1 are discussed in appendix B. In the same figures, we also plot — ^<I>„(x) using 
the scaling function ^n{x) obtained as described above. Then the solutions to (3.17) can 
be obtained as the x coordinates of the intersection points. 

For n = 1,4 there is only one intersection. The effect of the phase is to suppress the 
X ^ 1 region and to shift the corresponding scale to larger/smaller scales, respectively. For 
n = 2, 3 there are two intersections, which correspond to the local maxima of Pn{x)- The 
effect of the phase is to produce a double peak structure in pn{x), which corresponds to 
two dynamically generated length scales for the corresponding dimensions. In Table 2 we 
present the solutions of (3.17) for each n, where Xg and xi denotes the solutions in the 
X < 1 and X > 1 regions, respectively. 

^^We have checked that this is indeed the case in the region of x where 'w„{x) is close to 1. 
^*We will discuss in section 4.5 that in fact there should be more intersecting points that are not seen in 
these plots. The statements given below ignore this subtle issue. 
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Figure 2: Results for the r = 1 case. (Left) The function ■^logw„(x) is plotted for N = 
4,6,8 together with the scahng function obtained by large- extrapolations based on the 

asymptotic behavior (3.18). (Right) The function ■^fn\x) is plotted for AT = 16,32,64 together 
with — ^$„(a;). The x coordinate of the intersecting point gives a solution to (3.17). 
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Figure 3: Results for the r = 2 case. (Left) The function ■^logw„(x) is plotted for N = 
4,6,8 together with the scahng function obtained by large- extrapolations based on the 

asymptotic behavior (3.18). (Right) The function ■^fn\x) is plotted for AT = 16,32,64 together 
with — ^$„(a;). The x coordinate of the intersecting point gives a solution to (3.17). 
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Table 2: The first two columns show the solutions {xg and xi) to eq. (3.17) that correspond to the 
(local) maxima of Pn{x). The other columns show the results for $„(a;s), ^n{xi), S„ and A„ that 
appear in eq. (3.20). 

For n = 2, 3, we have to determine which of the two peaks of Pn{x) actually dominates 
in the A/" — > oo limit. For that purpose we calculate the difference of the free energy [24] 

An =^ -^|log/9„(xi) - logp„(xs)} = {$„(xi) - $„(Xs)} + S„ , (3.20) 

where S„ J\x | A/(o)(^)| . (3.21) 

If we find that A„ is positive (negative), the peak at x = xi (x = Xg) dominates in the 
large-A'' limit. Note that A„ depends on the values of ^nix) at Xg and x\, but there is 
no need to know the correct values of ^n{x) in the region where we have a severe sign 
problem. The integral S„ can be computed by fitting the data points for largest N to some 
known function with free parameters in the region Xg < x < X\, and then integrating that 
function over the region. In Table 2 we also present the results for $„(xs), $n(xi), H„ and 
A„. For n = 2 we obtain a clearly positive value for A„, which means that the peak at 
X = xi dominates. For n = 3 the obtained value for A„ is too small to reliably determine 
the dominant peak considering the expected systematic errors (See section 4.6.) in $„(xs) 
and ^n{xi)- 

The results presented above suggest that (A„) indeed deviate from 1 representing the 
strong effect of the phase, and that we are likely to get (Ai),(A2) > 1 > (A4), which is 
consistent with the GEM results. Thus our results support the speculation that the SSB of 
SO (4) symmetry occurs due to the effect of the phase. Moreover, the value of xi for n = 2 
and the value of Xg for n = 3 are reasonably close to the GEM results for (A2) and (A3), 
respectively, obtained with the S0(2) ansatz. (See Table 1.) 

On the other hand, our results for (Ai) and (A4) are not very close to the corresponding 
GEM results with the SO (2) ansatz. We also notice a puzzle. In usual Monte Carlo 
calculations, one obtains various observables using the same set of configurations. However, 
in the present method, one has to perform independent simulations to obtain (A„) for 
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different n. In fact the dominant configurations one obtains by simulating the constrained 
system (3.7) witli tlie parameter x fixed at tlie GEM results for (A^) can be quite different 
for different n. While this is not a problem per se at least theoretically, it suggests the 
existence of a practical problem, which is closely related to the observed discrepancies. 

For instance, the dominant configurations one obtains from the constrained system 
(3.7) with n = 1 and x > 1 have Ai(= x) > 1 > A2, A3, A4. Note that wi{x) is calculated as 
the VEV of e*'" in such an ensemble. When x is close to (Ai) obtained by the GEM with 
the S0(2) ansatz, the dominant contribution to wi{x) is expected to come from S0(2) 
symmetric configurations due to the enhancement by e*^. This clearly suggests that the 
overlap problem still remains. As a result, it is expected that wi{x) is underestimated in 
that region of x, which naturally explains the discrepancy for (Ai) observed above. Due 
to this remaining overlap problem, we need to be careful in interpreting the results of the 
single-observable analysis. We will come back to this point in section 4.5. 

4. Factorization method with multiple observables 

The important point of the method described in the previous section was to consider the 
distribution function (3.2) instead of trying to calculate the VEV directly by reweighting 
(2.16). Then the effect of the phase was factorized and represented nicely by the weight 
function Wn{x) defined by (3.6). This function suppresses the region near x = 1, which is 
preferred by the phase-quenched model, and enhances the regions x <^ 1 and/or a; S> 1. 
By solving the equation (3.17), we can obtain solutions quite different from x = 1. If we 
simulate the phase-quenched model, the observables An fluctuate around An = 1, and it 
would be difficult to sample configurations having A^ different from 1, which represents 
the overlap problem. By simulating the constrained system (3.7), we are able to estimate 
the effect of the phase in the region strongly suppressed in the phase-quenched model, and 
thus the overlap problem can be reduced. 

However, we also noticed that the overlap problem can still remain when one evaluates 
Wn{x) defined by eq. (3.6) since it is calculated as the VEV of e^^ in the constrained 
system (3.7). There might be some region of the configuration space which cannot be 
sampled efficiently by simulating the constrained system (3.7), and yet has important 
contribution to the VEV due to less fluctuation of F. In this section we discuss how we 
can solve this problem by generalizing the factorization method to multiple observables. 
Instead of considering the distribution function of A„ for each n separately, we consider the 
distribution function specifying all the An (n = 1, 2, 3, 4) at the same time. We discuss the 
possibility of including more observables in section 4.7. 

4.1 the generalized formulation 

We can actually generalize all the formulae in section 3.1 to the multi-observable case in a 
straightforward manner. The relevant functions for the single-observable case can then be 
written in terms of those for the multi-observable case. In this generalized formulation, one 
can clearly identify the overlap problem that still remains in the single-observable analysis. 
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Let us begin by defining the distribution functions for multiple observables as 

p{xi,X2, X3, X4) = ( JJ 5(aJit - Ajfc)^ , (4.1) 

k 

(xi, X2, X3, X4) = ( n ^i^k - Ait)) (4.2) 



for the full model and the phase-quenched model, respectively. By definition (2.8), these 
functions vanish unless xi > X2 > > X4. Applying the reweighting (2.16) to the right- 
hand side of (4.1), one finds that it factorizes as 

p{xi,X2,X3,X4) = p^^\xi,X2,X3,X4)w{xi,X2,X3,X4) , (4.3) 

where C is a normalization constant given by eq. (3.4). The correction factor w{xi,X2,xs, X4) 
is defined by 

w(^Xi, X2, X2, X4) = (e )xi,a;2,a;3,a;4 ~ {(^OsV^xi,X2,X3,X4, ) ('^•'^) 

where ( • )a;i,a;2,a;3,a;4 denotes a VEV with respect to the partition function 



7 



J dAe-^° Hdixk - Xk) . (4.5) 



Note that all the four observables A„ {n = 1,2,3,4) are constrained here in contrast to 
(3.7). 

We can rewrite the relevant functions in the single-observable case in terms of the 
functions defined above. For instance, we have 

Y[dXkp{xi,X2,X3,X4) , (4.6) 
fc^n 

Pn\xn) = J '[[dXkP^^\xi,X2,X3,X4) . (4.7) 

Using (3.3), we also obtain 

^ (^X )- - / Hfc^n dXk {Xl ,X2,Xs, Xj) w{xi , X2, X3, X4) 

Pn\xn) J'n.k^ndXkP^°Kxi,X2,X3,X4) 

where we have used (4.3), (4.6) and (4.7) in the second equality. This formula reveals the 
possibility that the overlap problem can still occur when one calculates the function Wn{x) 
defined by (3.6) as the VEV of e*^ in the system Zn,x- The region that gives dominant 
contribution in the numerator of (4.8) may not be sampled efficiently by simulating the 
system Zn^x- It is also clear that the overlap problem can be reduced further by constraining 
all the four observables A„ (n = 1, 2, 3, 4) at the same time. 

Now we have to maximize the distribution function p{xi,X2,xs, X4) with respect to xi, 
X2, X2 and X4. This leads to the coupled equations 

d d 

logp^'^\xi,X2,X3,X4) = -- — logi(;(a;i,X2,a;3,a;4) for n = 1, 2, 3, 4 . (4.9) 
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The function on the left-hand side and tt;(xi, X2, X3, 2:4) defined by (4.4) can be obtained 
similarly to the single-observable case described in section 3.2. In fact it is a formidable 
task to search for solutions in the 4d parameter space (xi, X2, ^3, X4) in full generality. 
However, given the insights we obtained from the analysis in the previous section, we can 
restrict ourselves to the region of (a;i, 0:2, a^a, 0:4) in which we expect to find a solution. In 
particular, it is expected that there exist solutions which satisfy 

(a) xi = X2 = X3 > 1 > a;4 , 

(b) Xi = X2 > 1 > Xs > X4 , 

(c) xi > 1 > a;2 > a;3 > X4 . (4.10) 

The solutions one obtains for the cases (a) and (b) correspond to the S0(3) and S0(2) 
symmetric vacua, respectively, and hence one can compare the results against those ob- 
tained by the GEM with the S0(3) ansatz and the S0(2) ansatz.-*^^ That would still require 
solving coupled equations for 2 variables in the case (a), and for 3 variables in the case 
(b). In what follows, we restrict ourselves to the r = 1 case for simplicity, and investigate 
whether the GEM results for (A„) with the S0(3) ansatz and the S0(2) ansatz are indeed 
solutions to (4.9) in the cases (a) and (b), respectively. We also discuss the calculation of 
free energy in section 4.6. 

4.2 analysis for the SO (3) symmetric vacuum 

Let us consider the case (a) in (4.10). We define the reduced functions 

Pso(3)(^'2/) =^ p''°\x,x,x,y) , (4.11) 
wso{3) {x, y) w{x, X, X, y) . (4.12) 
Prom the coupled equations (4.9), we obtain the counterpart of (3.17) as 

1 d 
J^fsoisui^'^y) = —Q^^so{3)ix,y) , (4.13) 

where C, = x,y, and we have defined 

42(3), c^^' = ^ log pfom y) ' (4-14) 

def 1 

^so(3)(a;,y) = ^l™^;^logw^so(3)(a;,y) ■ (4.15) 

The GEM results obtained for r = 1 with the S0(3) ansatz are (An)so(3) = 1-17 (^ = 2, 3) 
and (A4)go(3) = 0.5 (See Table 1.) We would like to check whether (x,y) = (1.17,0.5) 
indeed solves eq. (4.13). 



'^^The case (c) was not studied by GEM in ref. [35] since it is technically more difficult than the other 
cases due to the existence of more free parameters in the Gaussian action. However, our result A2 > 
discussed in section 3.4 strongly suggests that the solution for the case (c) has larger free energy than the 
solution for the case (b). 
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First let us consider the C, = x component of eq. (4.13). Let us set y to 0.5 and solve 
the equation for x to see if the solution agrees with the value 1.17. We can calculate 
w;so(3)(a;, 0.5) and /so(3) a:(^' ^■^) using the method described in section 3.2. In fig. 4 
(Left) we plot the function log wso(3) {x, 0.5) for N = 8, 12, 16. We make a large-iV 
extrapolation using the asymptotic behavior at large x to obtain $so(3) i^-i 0-5) represented 
by the solid fines. (See appendix D.l for the details.) Figure 4 (Right) shows that the 
solution lies at x = 1.151(2), which is consistent with the GEM value 1.17. 




Figure 4: (Left) The function logwgo(3)(2^, 0.5) is plotted against x for N = 8,12,16. The 
solid lines represent tlic function 4'so(3) 0-5) obtained by extrapolation to A'^ = oo as described 
in appendix D.l. (Right) The function 77t/so(3) a;(2;, 0.5) is plotted against a; for A'' = 16,32,64. 
The solid lines represent — ^$so(3)(a;, 0.5) obtained from the plot on the left. 

Next let us consider the C = y component of cq. (4.13). We set x to 1.17 and solve the 
equation for y. In order to obtain <I>go(3)(1.17, y), we have to make a large- extrapolation. 
For that purpose we first tried to use the asymptotic behavior of log ti;go(3) (1-17, y) at 
small y. However, it turned out that the finite- A?^ effects in log u)go(3) (1-17, y) at small y 
are much severer than in log u;so(3) (a^) 0.5) at large x investigated above.^^ We therefore 
use an "orthogonal" method here. 

Note that we were able to obtain the function $30(3) (x, y) at y = 0.5 from fig. 4 (Left). 
We can do the same thing for different y such as y = 0.45, 0.55. (See appendix D.l for the 
details.) In fig. 5 (Left) we plot $so(3)(l-17, y) for y = 0.45,0.50,0.55. By fitting the three 
points to a straight line, we obtain the derivative ^<I>go(3) (1.17, y) = —0.33(2) at y = 0.5. 
Figure 5 (Right) shows that the solution lies at y = 0.59(2), which is reasonably close to 
the GEM value 0.5 given the uncertainties involved in the analysis. 



^^The same trend can also be seen in the single-observable analysis relevant to the S0(3) ansatz. In fig. 2 
(Left), the finite- A/' elTects in logM;4(a:) at x < 1 are much severer than in logu'3(a:) at a; > 1. 
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Figure 5: (Left) The values of the ftinction $so(3)(l-17, y) at y = 0.45, 0.5, 0.55 are plotted against 
y. (Right) The function 773-/so(3) ^(l-l^, y) is plotted against y iov N = 16,32,64. The solid lines 
represent — ^$so(3)(l-17,y) at y = 0.5 obtained from the plot on the left. 

4.3 analysis for the SO(2) symmetric vacuum 

Let us consider the case (b) in (4.10). We define the reduced functions 

4o(2) y' ^) p'^^^ ^) ' (^-16) 

'Wso{2){x, y, z) w{x, X, y, z) . (4.17) 

Prom the coupled equations (4.9), we obtain the counterpart of (3.17) as 

1 d 

]v2 4o(2),c(^'y'^) = -Q^^&0{2){x,y,z) , (4.18) 

where ( = x,y,z, and we have defined 

/S0(2),c(^'2/'^) =^ ^log4o(2)(^'y'^) ' (4-19) 

dcf 1 

^SO{2){x, y,z) = l\m^j^logwso(2){x,y,z) . (4.20) 

The GEM results obtained for r = 1 with the S0(2) ansatz are (A„)so(2) = 1-4 {n = 
1,2), (A,3)so(2) = 0-7 and (A4)so(2) = 0-5 (See Table 1.) We would like to check whether 
{x,y,z) = (1.4,0.7,0.5) indeed solves eq. (4.18). 

First let us consider the = x component of cq. (4.18). We set y and z to 0.7 
and 0.5, respectively, and solve the equation for x. In fig. 6 (Left) we plot the function 
loguiso(2)(a^) 0-7, 0.5) for N = 8,12,16. We make a large- AT extrapolation using the 
asymptotic behavior at large x to obtain $go(2)(a^, 0.7, 0.5) represented by the solid line. 
(See appendix D.2 for the details.) Figure 6 (Right) shows that the solution lies at x = 
1.373(2), which is consistent with the GEM value 1.4. 

Next let us consider the Q = y component of eq. (4.18). We set x and z to 1.4 and 0.5, 
respectively, and solve the equation for y. In order to obtain $30(2) (1-4, y, 0.5), we have to 
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Figure 6: (Left) The function ■^logwso(2)(a;,0.7, 0.5) is plotted against a; for A'' = 8,12,16. 
The soUd line represents the function $so(2)(a^, 0.7, 0.5) obtained by extrapolation to A?" = oo as 
described in appendix D.2. (Right) The function -^f^^^^^^^lxjO.TjO.b) is plotted against x for 
N = 16,32,64. The solid lines represent — ^^'so(2)(a;, 0.7, 0.5) obtained from the plot on the left. 
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Figure 7: (Left) The values of the function $so(2)(l-4, y, 0.5) at y = 0.65,0.70,0.75 are plotted 
against y. (Right) The function ■/^/go(2) y, 0.5) is plotted against y for = 16,32,64. The 

solid lines represent — ^$so(2)(l-4,2/,0.5) at y = 0.7 obtained from the plot on the left. 



make a large- extrapolation. The region of y of the function logit;so(2)(l-4, y, 0.5) is 
restricted to 0.5 < y < 1.4, and we do not have any asymptotic behavior with respect to 
y that can be used for large-A^ extrapolation. Therefore we use the "orthogonal" method 
which was used in section 4.2 for studying the C, = y component of eq. (4.13). Since we were 
able to obtain ^*so(2)(a;) Vi z) for y = 0.7 and z = 0.5, we do the same thing for different y 
such asy = 0.65, 0.75. In fig. 7 (Left) we plot $so(2) (1-4, y, 0.5) for y = 0.65, 0.70, 0.75. By 
fitting the three points to a straight line, we obtain the derivative ^$go(2)(l-4, y, 0.5) = 
-0.243(4) at y = 0.7. Figure 7 (Right) shows that the solution hes at y = 0.649(4), which 
is consistent with the GEM value 0.7 given the uncertainties involved in the analysis. 

Finally let us consider the C = component of eq. (4.18). We set x and y to 1.4 and 
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Figure 8: (Left) The values of the function $so(2)(l-4,0.7, 2;) a.t z = 0.45,0.5,0.55 are plotted 
against z. (Right) The function 77?/go(2) 2(1-4,0.7,2) is plotted against z foi N = 16,32,64. The 
solid lines represent — ^$so(2)(l-4, 0.7, 2;) at 2; = 0.5 obtained from the plot on the left. 

0.7, respectively, and solve the equation for z. In order to obtain $30(2) (1-4, 0.7, z), we use 
the "orthogonal" method."^''' Since we were able to obtain $go(2) (a^j -2) for y = 0.7 and 
z = 0.5, we do the same thing for different z such as z = 0.45, 0.55. In fig. 8 (Left) we plot 
^>gO(2)(l-4, 0.7, z) at z = 0.45,0.50,0.55. By fitting the three points to a straight line, we 
obtain the derivative ^$so(2)(l-4, 0.7, z) = —0.357(4) at z = 0.5. Figure 8 (Right) shows 
that the solution lies &t z = 0.551(2), which is close to the GEM value 0.5. 

4.4 estimates on the systematic error 

One thing we notice from the results in sections 4.2 and 4.3 (Sec Table 3 for a summary.) 
is that the agreement with the GEM results is better for the large eigenvalues (A^) > 1 
than for the small eigenvalues (A„) < 1. Apart from the use of the "orthogonal" method 
in the latter case, we consider that the systematic error may be a possible reason. 

Let us first consider the SO (3) symmetric vacuum. The VEV (A„) is determined by 
solving (4.13). The left-hand side is determined accurately since there is no sign problem, 
and they behave as 

]^4S(3),x(^' 0-5) - -«(^ - 1-17) -b, a = 1.43(2) , b = 0.174(1) , (4.21) 

]^/sS(3),y^l-l^'^) ^ -a(i/-0.5)+6, a = 1.70(4) ,6 = 0.452(1) , (4.22) 

near the solution as one can obtain from the right panel of figs. 4 and 5. Let us denote 
the right-hand side of the ( component of (4.13) by i^so(3),C' ^^'^ assume that it has an 
error A£)go(3),ij, where ( = x,y. This error includes the systematic error due to deviation 
from the asymptotic behavior as well as the one due to the large-AT extrapolation. Then 

^'^We first tried to use the asymptotic behavior of log wso(2)(l-4, 0.7, z) at small z. However, it 
turned out that the finite- Af effects in log wso(2) (1-4, 0.7, z) at small z are much severer than in 
loguiso{2)(ic, 0.7, 0.5) at large x investigated above. 
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the error in i?so(3),a; propagates to the error of the solution x as a Ax = ADgo(3),a;, from 
which we obtain 

Ax ^ |^SO(3),xl A.DsO(3),a; ^ ^ ^DsO{3),x ^^23) 
X ax |-DsO(3),a;| |-DsO(3),xl 

On the other hand, the error in -Dgo(3),y propagates to the error of the solution y as 
a Ay = AI?go(3),y, from which we obtain 

Ay _ \DsO{3),y\ ^Dso{3),y ^ ^ ^ ^DsO{3),y ^4) 



y ay |-Dso(3),2/| \Dso(3),y\ 

Thus we find that the coefficient is four times bigger for y than for x. A similar analysis 
can be made for the SO (2) symmetric vacuum, and we obtain 

Ax^ 0.24 AZ^so(2),. ^ p 2^ AI^so(2),. 

X 0.82 • 1.4 \Dso^2),x\ ■ I^SO(2),xl ' ^ ' ' 

Ay^ 0.25 Aj^so(2),. ^ AZ^so(2).. 

y 0.9-0.7 |I)so(2),,l ■ l^so(2),,l ' ^ ' ^ 

A. ^ 0-35 AZ^so(2).. ^ Q_43 AZ^so(2).. ^ 

z 1.6 • 0.5 |-Dso(2),zl l^so(2),2l 

where £^so(2),c AL>so(2),^ are defined analogously to the S0(3) case. Thus we find 
that the relative errors for the large eigenvalues tend to be smaller than those for the small 
eigenvalues. 

4.5 comparison with the single-observable analysis 

In this section wc compare the results obtained above by the multi-observable analysis 
with those obtained by the single-observable analysis in section 3.4. For that, we first need 
to reconsider the interpretation of the latter taking into account the artifacts due to the 
remaining overlap problem. 

Let us recall that the generalized distribution function /9(xi, X2, X3, X4) is expected to 
have three local maxima corresponding to each of the three cases (4.10). Therefore, the 
distribution function pn{x) in the single-observable analysis is also expected to have three 
local maxima due to the relation (4.6). The fact that we observe only two for n = 2, 3 and 
one for n = 1,4 as we described in section 3.4 is due to the remaining overlap problem. 

We reconsider the plots in the right column of fig. 2 from this point of view. For 
instance, the intersection of the two curves shown in the uppermost panel can be naturally 
identified as the one corresponding to the case (c). As x decreases from the intersecting 
point, the true curve y = "^^1(2^) is expected to have four more intersections with the 
curve y = -^f\^\x) in the x > 1 region. Two of them represent the local maxima of 
the distribution function pi(x) corresponding to the cases (a) and (b) in (4.10). Such 
a structure is invisible in the results of the single-observable analysis due to the overlap 
problem. From this point of view, it is misleading to compare the solution x\ for the n = 1 
case with the (Ai) for the 80(2) symmetric vacuum. Rather we should consider it as an 
estimate for (Ai) in the vacuum corresponding to the case (c) in (4.10). This provides a 
prediction for the calculation using GEM, which is not done yet. (See footnote 15.) 
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Similar arguments apply to the n = 2, 3, 4 cases. The solutions Xg and xi for the n = 2 
case should be interpreted as estimates on (A2) in the vacuum corresponding to the case 
(c) in (4.10) and in the S0(2) symmetric vacuum, respectively. The solutions Xg and x\ 
for the n = 3 case should be interpreted as estimates on (A3) for the S0(2) and S0(3) 
symmetric vacua, respectively. The solution Xg for the n = 4 case should be interpreted as 
an estimate on (A4) for the SO (3) symmetric vacuum. 

In Table 3 we summarize our results for (A„) corresponding to the SO (2) and SO (3) 
symmetric vacua. We also show the results of the single-observable analysis with the new 
interpretation except for an estimate on (A4) for the SO (2) symmetric vacuum, which is 
not available. We do find that the results of the multi-observable analysis have better 
agreement with the GEM values. 



ansatz 


S0(3) 


S0(2) 


method 


single-obs. 


multi-obs. 


GEM 


single-obs. 


multi-obs. 


GEM 


(Ai) 






1.17 






1.4 


(A2) 






1.17 


1.317(1) 


1.373(2) 


1.4 


(A3) 


1.129(1) 


1.151(2) 


1.17 


0.62(2) 


0.649(4) 


0.7 


(A4) 


0.71(5) 


0.59(2) 


0.5 


not available 


0.551(2) 


0.5 



Table 3: The results for the normalized eigenvalues (A„) for r = 1 obtained by the factorization 
method with the single-observable and multi-observable analyses for the SO (3) and SO (2) symmetric 
vacua. The dash implies that the result should be the same as the one below in the same column 
due to the imposed symmetry. We also show the GEM results obtained at = 00 in ref. [35]. 



4.6 calculation of the free energy 

In sections 4.2 and 4.3, we have checked that the GEM results with the S0(3) ansatz and 
the SO (2) ansatz can indeed be obtained as solutions to (4.9). In order to determine which 

solution dominates the path integral, we need to compare the free energy. This was done in 
the GEM and the result is shown in Table 1, from which we find that the SO (2) symmetric 
solution dominates in the large- limit. Here we try to obtain the difference of free energy 
for the two solutions by Monte Carlo simulation. The factorization method has the nice 
feature that the difference of free energy is decomposed into the term coming from the 
phase-quenched model and the term due to the effect of the phase. The former can be 
calculated accurately, and the latter is the main source of the error. 

What we should calculate is the difference of log p{xi, X2, X3, X4) between the two so- 
lutions. Let us denote the solution for cases (a) and (b) as Xa = {X' ,X' ,X' ,Y') and 
Xh = {X, X, Y, Z) , respectively. From Table 1 we find for r = 1 that 

(a) X' ~ 1.17 ,y' ~ 0.5 , (4.28) 
(6) X ~ 1.4 ,F ~ 0.7 ,Z ~ 0.5 . (4.29) 
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Using the factorization property (4.3), we obtain 

^ =^ ]^{log/9(Xa) -logp(fb)} 

= [<^soi3)iX',Y') -'^SOi2)iX, Y,Z)]+E, (4.30) 
where ^"^^ J dxj^-^logp^°\xi,X2,X3,X4) ■ (4.31) 

The first term in (4.30) can be estimated as (See appendix D.) 

$go(3) (1.17, 0.5) = -0.160(3) , (4.32) 
$go(2)(1.4,0.7,0.5) = -0.155(1) . (4.33) 

In order to calculate the second term H in (4.30), wc first define a path in the (zi, X2, X3, X4)- 
space connecting the two solutions Xa and Xf,, obtain the gradient of \og p^^\xi, X2, X3, X4) 
and integrate it along the path. As a path connecting the two solutions, we consider 



Ci 
C2 
C3 



{X,X,Y,Z)^{X,X,X,Z) , 
{X, X, X, Z) ^ {X, X, X, ¥') , 
{X,X,X,Y')^{X',X',X',Y') 



The paths Ci and C3 are investigated already in sections 4.3 and 4.2, respectively, while 
the path C2 should be investigated newly. In the r = 1 case, we have Y Z' accidentally, 
and therefore the study of the path C2 can be totally omitted. The second term in (4.30) 
can then be evaluated as 



1 AO) 



(1.4, y, 0.5) dy-J^ '^^ ]^/sS(3),.(^' 0-5) dx 



^~ J0.7 iV2^0(2),?; 

= (-0.014) - (-0.079) = 0.065 , (4.34) 

where the integrals are calculated by fitting the data points for largest N in figs. 4 and 7 
(Right) to some known functions with free parameters in the integration domain. The errors 
are negligible compared to those in the first term of (4.30). Thus we obtain A = 0.060(4), 
which should be compared with the value A c± —0.3 predicted by GEM as one can see from 
Table 1. 

According to the interpretation of the single-observable analysis given in section 4.5, 
what we have done above may be regarded as a refined version of the calculation of A3 in 
section 3.4, where we obtained A3 = 0.11(4) for r = 1. 

Let us recall that the calculation of (4.34) does not have much ambiguity. Note also 
that S > implies that the S0(3) symmetric configuration {X' , X' , X' ,Y') has lower free 
energy in the phase- quenched model than the S0(2) symmetric configuration {X, X, Y, Z). 
This is reasonable considering that the former configuration is closer to the dominant SO (4) 
symmetric configuration A„ = 1 in the phase-quenched model. 

On the other hand, the estimates (4.32) and (4.33) may be subject to systematic 
errors due to the deviation of the functions $so(3)(2^) 0-5) and ^so{2){xj 0.7, 0.5) from their 
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asymptotic behaviors and due also to large-A^ extrapolations. Let us recall that the error 
propagation in solving the = x component of eq. (4.13) is given by (4.23), while the error 
in (4.32) can be roughly estimated to be 

A^so(3)(l-17,0.5) ADgo(3),, 

|$so(3)(l-17, 0.5)1 ~ \Dsoi3),x\ ' ^ ' ' 

Therefore, the relative error in $go(3) (1-17, 0.5) is ten times larger than that in (A3) in the 
S0(3) symmetric vacuum. The crucial point here is that the function ^so(3)(^) 0.5) changes 
very rapidly in the region where the solution to (4.13) lies. Therefore, a small systematic 
error in its estimation affects the value at the solution much more than the solution itself. 
A similar argument applies also to the S0(2) symmetric vacuum. It is therefore reasonable 
to conclude that the first term in (4.30) is not obtained with the accuracy that is sufficient 
for the determination of the sign of A. 

While this sounds a bit disappointing, let us recall that the S0(2) and S0(3) sym- 
metric vacua become degenerate as the parameter r = Nf/N of the model goes to zero. 
For instance, it would have been much easier to compare the free energy for the solutions 
corresponding to the cases (b) and (c) in (4.10) as the calculation of A2 in section 3.4 sug- 
gests. Clearly the situation would be model dependent, and we consider that the difficulty 
in the calculation of the free energy is not a generic feature of the factorization method. 

4.7 including more observables 

As we have seen above, the overlap problem can be reduced drastically by constraining all 
the An (n = 1, 2, 3, 4) at the same time instead of constraining just one of them. A natural 
question that arises then is whether there is no more overlap problem we have to worry 
about so that we actually do not need to include more observables in our analysis. Note 
that this is relevant not only to the calculation of the VEV of observables other than A„, 
but also to the calculation of (A^) and the free energy that we have discussed above. 

Let us see directly a possible overlap problem associated with a general observable O 
in a simulation with constraints on A„ (n = 1, 2, 3, 4). For that we rewrite the VEV (O) as 



» 4 4 

(O) = / n dxk (p n 5{xk - \k)) (4.36) 
k=i k=i 

4 (oeA 



fe=i \c' , 

\ / X\,X2,X3,X4 

Assuming the usual equivalence between the canonical ensemble and the microcanonical 
ensemble, we may consider that the integral over Xk in (4.37) is dominated by (X, X, Y, Z) 
giving the absolute maximum of /9(xi, X2, X3, X4). This leads to 

(p) ~ ' ' x.x.y.z _ (4_3g) 
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Thus we can calculate the VEV by using the standard reweighting formula within the 
"microcanonical ensemble" characterized by {X, X,Y, Z). For O = \n, the VEV will be 
given by X, X, Y, Z for n = 1, 2, 3, 4, respectively, as it should. 

For a general operator O, however, we can still have an overlap problem due to the 
same reason that we have one when we apply the reweighting method to the original model. 
If it turns out that {0)x,x,Y,z is actually close to (4.38), there is no more overlap problem 
as far as the operator O is concerned. Note that this requires that 

Oc'^) r^(o) (e^^) , (4.39) 

/ X,X,Y,Z \ I X,X,Y,Z \ I X,X,Y,Z 

as we can sec from (4.38). Eq. (4.39) only implies that the correlation of the operator 

O and the phase factor e*^ within the "microcanonical ensemble" is small. If the same 

statement holds for arbitrary observable, we may say that the sign problem is practically 

solved even if ( e*^ ) is not close to one at all. 

\ I X,X,Y,Z 

Thus we find that the success of the factorization method in a general model relies 
on whether one can find the minimal set of observables that is needed to "solve" the sign 
problem in the above sense [36]. We consider that the matrix model we are studying 
provides an explicit example in which this can be done. 

In order to provide some evidence for this statement, let us consider an observable 

= -^J2^v[A„A,]\ (4.40) 

and the normalized one O = O /{0)q. When we constrain O to a small value, the dominant 
configurations have the property [A^,74i/] « 0, meaning that are simultaneously diag- 
onalizable, e.g. as = diag(a^^\ • • • ,aj^^). For such configurations, the determinant 
becomes 

N 

detP = n{E(«?)'}^0- (4-41) 

Therefore, the observable (4.40) is considered as a "dangerous" one, which can potentially 
have strong correlation with the phase factor. 

Clearly the formula (4.38) is not useful for actually investigating the overlap problem 

associated with the operator O. For that purpose, we can use the factorization method 
including not only A„ (n = 1, 2, 3, 4) but also O. Let us define the corresponding functions 
p^'^) and w with five arguments, and also define the reduced functions 

pg)(x) = p(0)(X,X,F,Z,x) , (4.42) 
wo{x) = w{X, X, Y, Z, x) . (4.43) 

The VEV (O) can then be obtained by solving 

W2foi^) = -^^oix), (4.44) 
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Figure 9: (Left) The ftinction -^logwoix) is plotted against a; for A'' = 6,8,12. We also plot 
the function $0(2;) obtained by extrapolation to iV = 00 as described in appendix D.2. The 
two solid lines represent the margin of error. (Right) The function ^ log (a;) is plotted for 
N = 8, 16, 32. We also plot —^^"0(2;) obtained from the plot on the left. 

where fQ\x) and ^o{x) are defined as 

f§\x) = ^\ogp%\x), (4.45) 
$o(x) = lim -^logwoix) . (4.46) 

AT— >-oo iV 

Prom fig. 9 (Left) we find that -^logwoix) approaches zero for x — > as expected 
from (4.41). We make a large-AT extrapolation as described in appendix D.2 to obtain 
^o{x), which is also shown in the same figure by the two solid lines showing the margin 
of error. From fig. 9 (Right) we find that the effect of the phase is to shift the estimate 
of (O) by Ax = 0.07(3). On the other hand, the standard deviation of the distribution 
Pq{x) is estimated as cr ~ 0.7 /N from the slope of the function plotted in fig. 9 (Right) 
around x ~ 0.92. This means that the deviation Ax is < 2a for N < 16. Thus, the 
remaining overlap problem associated with this observable (4.40) is practically small. This 
is consistent with the fact that we were able to reproduce the GEM result by constraining 
only the four observables (n = 1, 2, 3, 4). 



5. Summciry and discussions 

In this paper we tested the proposed scenario for dynamical compactification of space- 
time in the IIB matrix model, in which the phase of the fermion determinant induces the 
spontaneous breaking of rotational symmetry. We have shown that this can indeed occur 
in a simplified model [34] by performing Monte Carlo simulations. The model exhibits SSB 
of SO (4) rotational symmetry in agreement with the results of GEM [35], and the phase 
of the fermion determinant plays a crucial role in the mechanism of SSB. 

First, in the absence of the phase, we have confirmed that the model has no SSB 
by calculating the VEV of the eigenvalues of T^^, which is analogous to the moment of 
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inertia tensor. The effect of the phase has been studied by using the factorization method 
originally proposed in ref. [24]. From the single-observable analysis, we find that the phase 
fluctuations strongly suppress the region of the configuration space favored in the phase 
quenched model and result in different length scales for each dimension of space-time. While 
our results are in partial agreement with the GEM results in ref. [35] , we also observe some 
discrepancies and puzzles, which wc attributed to the remaining overlap problem. This 
motivated us to generalize the method to multiple observables. In particular, we find 
that controlling the four eigenvalues of the "moment of inertia tensor" is good enough to 
reduce the overlap problem to a sufficient level. Restricting ourselves to the SO (2) and 
S0(3) symmetric vacua, we were able to reproduce the GEM results for the VEV of the 
eigenvalues consistently. 

Our results are an encouragement for pursuing similar studies on the IIB matrix model. 
Although it is computationally more demanding, the supersymmctry of the IIB matrix 
model could make the SSB of rotational symmetry easier to see. In particular, we have 
observed in ref. [24] that the ■^fn\x) that appears in (3.17) is actually very close to zero, 
which can be understood as a result of cancellations by fermionic and bosonic contributions 
to the interactions among space-time points. This implies that the solutions to cq. (3.17) 
appear in the region where the sign problem is not severe. We are currently studying the 
6 dimensional version of the IIB matrix model by Monte Carlo simulation [37] and trying 
to compare the results against the predictions obtained by the GEM recently [38] . 

We would like to emphasize that the model studied in this paper has a severe sign 
problem despite its simpleness. It is encouraging that such a system can be studied by 
Monte Carlo simulation using the factorization method, which gives us a lot of useful 
insights into the effect of the phase. The method itself is quite general and it can be 
applied to any interesting system which suffers from the sign problem. The crucial step 
is to find an appropriate set of observables that one can control in order to determine 
and sample efficiently the region of the configuration space favored by the relatively small 
fluctuation of the phase. We expect that this is possible in many interesting systems. 
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A. Details of the Monte Carlo simulation 

In this section wc present the details of the algorithm used for our Monte Carlo simulation. 
It is essentially the hybrid Monte Carlo (HMC) algorithm, which has been applied in similar 
models in refs. [21,22,24,25,39]. The computational effort grows as 0(A^^) in the present 
model. 
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We first discuss tlie HMC algoritlim for tlie phase-quenclied model (2.12). We introduce 
auxiliary bosonic Hermitian matrices and consider the action 

Shmc[P, A] = ^ tr(P^)2 + So[A] . (A.l) 

The original model is obtained by integrating out P^. We regard the action iShmc [-P; ^] 
as the Hamiltonian of a classical system described by ^^(r) and its conjugate momentum 
Ph{t), where r denotes the fictitious time of the classical system. The Hamiltonian equation 
of motion is obtained as 



diA^ _ dSuMC _ ^ r, ^ /A 9^ 

dr d{Py,)ki 

d{P^)ki ^ _dSBMC_ 
dr d{An)M 

,..3) 

where Tr denotes a trace with respect to the 2A?^-dimensional index. The indices k, I run 
over k,l = 1,2, ■■ • ,N. The updating procedure consists of (i) refreshing the momentum 
variables by Gaussian random numbers, which obey the distribution e~'^HMc ^^d (ii) 
solving the Hamiltonian equation of motion for a fixed time interval r. In actual calculations 
we have to discretize the Hamiltonian equation (A. 2) and (A.3). The reversibility of the 
time evolution is preserved by using the so-called leap-frog discretization, which gives 



(AW)« = (4°))^, + At {Pj^'^^k , (A.4) 
(4''^'^)^^ = (4% + Ar (p(''+i/2))zfc , (A.5) 

where r=l,2, ••• ,1^ — 1 and T = v Ar, and we have introduced the short-hand notation 
pj[^ = P^(rAr) and 4^ = ^/i('^Ar). The conservation of the Hamiltonian (5hmc) is 
violated by the discretization. The detailed balance can be preserved, however, by adding 
a Metropolis accept/reject procedure at the end of each trajectory. Namely we accept 
the trial configuration with the probability min(l, e~^'^HMc)^ where AS'hmc represents the 
difference of S'hmc between the trial and original configurations. The step size Ar for the 
time evolution should be small enough to keep the acceptance rate reasonably high. 

Next wc discuss how to simulate the constrained system defined by the partition func- 
tion Zxi,x2,x3,x4 ill eq. (4.5). As we discussed in section 3.2 for the single-observable analysis. 
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we approximate the delta function by the Gaussian function as 

(so + ^VniXn) 



Zxi,x2,x3,x4,,v — j dA exp 



n=l 



(A.7) 



Vn{K) = \ln{\n-inf . (A.8) 

The system (3.10) used for the single-observable analysis in section 3 can be obtained by 
setting 7fc = for 7^ n. 

Taking into account the potential term Vn{\n) in (A.7), we need to subtract the term 
^diA^)ki right-hand side of eq. (A. 3). This term can be calculated explicitly as 

follows. Let us note first that the eigenvalues A„ satisfy 

4 

Y^T.pvf^ = Xnv'^\ (A.9) 

p=l 

where v^^^ is the eigenvector of the 4x4 matrix T^^ normalized as X]/i=i '^'m^^^A'"^ ~ 
summation over the index n.) Taking the derivative of eq. (A.9) with respect to (A^)jk/, 
we obtain 

^ \d{An)ki ^ d{An)ki ) d{Ai^)ki d{A^)ki 

Multiplying both sides of eq. (A.IO) by and taking a sum over i/, we obtain 

where the second terms of each side of eq. (A.IO) cancel. Therefore we obtain 

||^ = 7n(A„-C„)^ = ^(A.-Ui:<M")(^^k ■ (A.12) 

Let us comment on actual values we have chosen for the parameters in the potential 
(A.8). For the single-observable analysis in section 3, we typically^^ used 7„ = 10^. For 
the multi-observable analysis in section 4, we studied the cases (a) and (b) in (4.10) cor- 
responding to the SO (3) and SO (2) symmetric vacua, respectively. In the case (a), we set 
7i = 72 = since Ai and A2 are automatically close to A3. When wc fix A4 to the value 
obtained by GEM and vary A3 as we do in section 4.2, we set 74 much larger than 73 
(typically 74 = 10^ and 73 = 10^), and set ^4 to the GEM value for (A4) obtained with the 
SO (3) ansatz. In the case (b), we set 71 = since Ai is automatically close to A2. When 
we fix A3 and A4 to the values obtained by GEM and vary A2 as we do in section 4.3, we 
set 73 and 74 much larger than 72 (typically 73 = 74 = 10^ and 72 = 10^), and set ^3 and 
^4 to the GEM values for (A3) and (A4), respectively, obtained with the S0(2) ansatz. 



^**In fact we have used 7„ in the range 10^ — 10^. We have checked that our results are independent of 
7„ within a wide range. 
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B. Asymptotic behaviors of the functions fn\x) 

In this section we discuss the asymptotic behaviors of the functions fn\x) that appears 
in (3.17). While these behaviors are not exphcitly used in our analysis, they help us 
understand the dynamical properties of the eigenvalues A„ in the matrix model. In fact 
the functions fn\x) in the IIB matrix model behaves quite differently as we mention in 
section 5, and wc hope that our results given below would be useful for comparison. 



Let us recall that the function fn \x) is defined by eq. (3.9) in terms of pn' {x), which 
is the distribution function of A„ in the phase-quenched model. It is therefore expected to 
have simple scaling behavior when x <C 1 or x S> 1. At small x we expect 



(0), 



n) + rSni > - + an 

X 



(B.l) 



This behavior is due to the phase space suppression since (5 — n) directions shrink as we 
decrease x. Considering that the shrunken directions have the extent proportional to ^/x, 
wc obtain pn\x) ~ (^/x)^^^^'^^^'^ . The n = 1 case differs since all the eigenvalues of 
collapse to at a; <C 1 and the suppression factor comes also from the fermion determinant, 
which is a homogeneous polynomial of of degree 2rN^. This gives an extra suppression 
factor of {■s/x)'^'^^^ . Prom the definition (3.9), we obtain the leading behavior of (B.l). 
Assuming that we have some analytic function of x multiplied to the leading power-law 
behavior of pn\x), we obtain some constant a„ in (B.l) as the subleading term. 
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Figure 10: The function -^fn [x) at iV = 64 is plotted against x for r = 1 (Left) and r = 2 
(Right). The straight lines are fits to the behavior (B.l). 

At large x we expect 

^4'^)W--^n(A„)o + g+r)l. (B.2) 

As we increase x, the eigenvalues Ai,--- , A„ are forced to be large and the (n + l)-th 
through the 4-th directions remain relatively small. In that case it is the bosonic action 
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Figure 11: The function ■^fn\x) at A'' = 64 is plotted against 1/a; for r = 1 (Left) and r 
(Right). The straight hnes represent the behavior (B.2). 



= 2 



that dominates the suppression of those configurations and we expect that 



PnHx) - exp { - 



iiVtr (A^)2)^ J ^ exp |^-iiv2f;(A,)„,v^ 



(B.3) 



Since {Xk)n,v ~ x{\n)o at large x foi k = 1, ■ ■ ■ ,n, we obtain the first term of (B.2), which 
becomes —^n{l + |) at large N. We also expect a power-law correction (v^)("'+^'')^^ to 
pn^ (x) that comes from the measure and the fermion determinant as we discussed in the 
case of a; ^ 1. Thus we obtain the second term in (B.2) as the subleading term. Figures 
10 and 11 confirm the above asymptotic behaviors at small x and large x, respectively. 



C. Large- A'^ extrapolations in the single-observable analysis 

In this section we explain the large- AT extrapolations made in the single-observable analysis 
based on the asymptotic behavior (3.18). 

In figs. 12 and 13 we show the log- 
log plots of log Wn{x) for r = 1 and 
r = 2, respectively. We fit the data to the 
behavior (3.18) at x ^ 1 and at x S> 1, 
and extract the coefficients Cn and cZ„ in 
(3.18) for each N. The results are plot- 
ted against jf in figs. 14 and 15 for r = 1 
and r = 2, respectively. Assuming that 
the finite- AT effects are 0(1/A^), we make 
large- AT extrapolations, which give the val- 
ues shown in Table 4. The corresponding 
functions (3.18), which wc define as the scaling functions ^nix) in each region, are also 
plotted in figs. 12 and 13. (The two solid lines show the margin of error.) The scaling 
functions thus obtained are plotted in the left column of figs. 2 and 3. 





r 


= 1 


r = 


= 2 


n 




dn 






1 




2.26(6) 




2.02(5) 


2 


2.3(1) 


0.378(2) 


2.407(2) 


0.53(2) 


3 


0.52(4) 


0.20(3) 


0.96(6) 


0.353(3) 


4 


0.26(4) 




0.46(8) 





Table 4: Large-A/' extrapolated values of the co- 
efBcients c„ and dn in the asymptotic behaviors 
(3.18). 
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Figure 12: The log-log plot of — -^logu>„(x) for r = 1. The straight solid lines represent the 
power-law behavior (3.18) with the coefficients presented in Table 4, which are obtained by the 
large- extrapolation from the N = 4,6,8 data as described in fig. f4. A clear trend towards 
large- A'' scaling is observed, in particular, for the n = 2, 3 cases. 
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Figure 13: The log-log plot of — -^logu>„(x) for r = 2. The straight solid lines represent the 
power-law behavior (3.18) with the coefficients presented in Table 4, which are obtained by the 
large- extrapolation from the N = 4,6,8 data as described in fig. 15. A clear trend towards 
large- A'' scaling is observed, in particular, for the n = 2, 3 cases. 
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Figure 14: The coefficients c„ and c?„ in the asymptotic formula (3.18) for r = 1 extracted from 
fig. 12 are plotted against jf. The data points are fitted to a straight line. The large- iV extrapolated 
values obtained in this way are presented in Table 4. 
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Figure 15: The coefficients c„ and c?„ in the asymptotic formula (3.18) for r = 2 extracted from 
fig. 13 are plotted against jf. The data points are fitted to a straight line. The large- iV extrapolated 
values obtained in this way are presented in Table 4. 
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D. Large- A*" extrapolations in the multi-observable analysis 

In this section we summarize the laxge-iV extrapolations we made in the multi-observable 
case. We describe them for the S0(3) and S0(2) symmetric vacua separately. 

D.l extrapolations for the SO(3) symmetric vacuum 

Similarly to the asymptotic behavior (3.18) used in the single-observable analysis, it is 
expected that 

log u;so(3) (x, 0.5) - -dix-^ + d2X-^/^ (D.l) 

at X S> 1. Here we have added the subleading term motivated from the fact that oc 
l/-\/x as described below (3.18). 

In fig. 16 (Top), we plot log u'go(3)(x, 0.5) against which confirms (D.l) includ- 
ing the subleading term. We plot the coefficients di and d2 against in fig. 16 (Bottom), 
which shows that the finite N effects are of the order of 1/A^. Based on this observation, 
we obtain the large- extrapolated values di = 0.1572(5) and d2 = —0.032(4). The scaling 




Figure 16: (Top) The function log wso(3)(a;, 0.5) is plotted against for N = 8,12,16. The 

straight lines represent the fits to the asymptotic behavior (D.l). (Bottom) The coefficients di and 
d2 in the asymptotic behavior (D.l) extracted from the top figure are plotted against j^. The data 
points can be fitted nicely to a straight line. 
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function (4.15) obtained in this way is plotted in fig. 4 (Left). By plugging x = 1.17 into 
(D.l), we obtain $go(3)(1.17, 0.5) = -0.160(3) as in eq. (4.32). 

We redo the calculation for y = 0.45, 0.55 and obtain the large-iV extrapolated values 

di = 0.1433(2) , d2 = -0.023(1) for y = 0.45 , (D.2) 
di = 0.1749(2) , d2 = -0.029(3) for y = 0.55 , (D.3) 

from which we obtain $go(3)(1.17, y) for y = 0.45 and 0.55. These results, together with 
the value at y = 0.5, are plotted in fig. 5 (Left). 

D.2 extrapolations for the SO(2) symmetric vacuum 

Similarly to (D.l) for the S0(3) symmetric vacuum, it is expected that 

log u;so(2) {x, 0.7, 0.5) ~ -dix'^ + d2X-^/^ (D.4) 

at X ^ 1. In fig. 17 (Top), we plot log^i'so(2) (a^j 0.7, 0.5) against which confirms 
(D.4) including the subleading term. We extract the coefficients di and d2 for N = 8, 12, 16, 
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Figure 17: (Top) The function logwso(2)(a;,0.7,0.5) is plotted against fov N = 8,12,16. 
The straight lines represent the fits to the asymptotic behavior (D.4). (Bottom) The coefficients di 
and d2 in the asymptotic behavior (D.4) extracted from the top figure are plotted against j^. The 
data points can be fitted nicely to a straight line. 
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and make a large- extrapolation for di and d2 as in fig. 17 (Bottom), from which we get 
the large-A*" extrapolated values di = 0.322(2) and d2 = 0.021(1). The scaling function 
(4.20) obtained in this way is plotted in fig. 6 (Left). By plugging x = 1.4 into (D.4), we 
obtain $so(2)(l-4,0.7,0.5) = -0.155(1) as in eq. (4.33). 

We redo the calculation for y = 0.65,0.75, and obtain the large- AT extrapolated values 

di = 0.295(3) , d2 = 0.018(1) for y = 0.65 , (D.5) 
di = 0.344(3) , d2 = 0.019(3) for y = 0.75 , (D.6) 

from which we obtain ^>so(2) (1-4, y, 0.5) at y = 0.65, 0.75. These results, together with the 
value at y = 0.7, are plotted in fig. 7 (Left). 

We redo the calculation for z = 0.45, 0.55, and obtain the large- A*" extrapolated values 

di = 0.292(5) , d2 = 0.026(2) for z = 0.45 , (D.7) 
di = 0.352(2) , d2 = 0.015(3) for z = 0.55 , (D.8) 

from which we obtain <I>so(2)(l-4, 0.7, z) at z = 0.45, 0.55. These results, together with the 
value at 2; = 0.5, are plotted in fig. 8 (Left). 



-0.05 

S -0.1 

o 

I, -0.15 
o 

C -0.2 



N=6 - 




J 


N=8 






N=12 - 





































0.2 0.4 0.6 0.8 1 




0.05 0.1 0.15 0.2 0.05 0.1 0.15 0.2 

1/N 1/N 



Figure 18: (Top) The function -^^tj^ logw;o(a;) is plotted against ^Jx for = 8, 12, 16. (Bottom) 
The coefRcients d\ and di in the asymptotic behavior (D.9) are plotted against j^. The data points 
can be fitted nicely to a straight line. 
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Finally let us describe the large-A^ extrapolations used in section 4.7 in the analysis 
for the observable (4.40). First the asymptotic behavior of log wo{x) is expected to be 

^ log wo (x) = -dix'^ + d2X^/^ . (D.9) 

The leading term is obtained as follows. Let us consider a small perturbation around a 
diagonal configuration. The phase of the fermion determinant appears at the second order 
of this perturbation. On the other hand, the observable O becomes nonzero also at the 
second order. Therefore the distribution of the phase is expected to have a width a oa x. 
Applying the formula (3.19), we obtain the leading term. The power of the subleading term 
can be deduced by expanding the fermion determinant around a diagonal configuration. 

In fig. 18 (Top) wc plot ^^^^ log ti;c)(x) against ^/x. The data points can be fitted 
nicely by straight lines, which confirms the asymptotic behavior (D.9). We extract the 
coefficients di and d2 for = 6, 8, 12, and make a large-A extrapolation for di and d2 as 
in fig. 18 (Bottom). We obtain di = 0.38(2) and (I2 = 0.26(1). The scaling function (4.46) 
obtained in this way is plotted in fig. 9 (Left). 
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